# Code to make Figures for Mattes and Weeks IO "Apology Diplomacy: The International Image Effects of Interstate Apologies"

################# Create Figure 1 ################### 

# Clear environment
rm(list = ls())

# Load necessary libraries
library(dplyr)
library(ggplot2)
library(broom)
library(haven)
library(cowplot)
library(devEMF)
library(png)
library(grid)
library(labelled)

setwd("~/Dropbox/HawksDoves/Apologies/Apologies--Image/Submission/IO/Resubmit/Conditional Accept/Replication Data/R code")

# Load the data
data_full <- read_dta("Study1_prepped.dta")

# Exclude "long" apologies
fig1data <- data_full %>% filter(apol_type != "long")

# Function to run regression and tidy results
run_regression <- function(dv, data) {
  model <- lm(as.formula(paste(dv, "~ autnone + autshort + demnone + demshort - 1")), data = fig1data)
  tidy_model <- tidy(model, conf.int = TRUE)
  tidy_model$term <- factor(tidy_model$term, levels = c("autnone", "autshort", "demnone", "demshort"))
  list(model = model, tidy_model = tidy_model)
}

# Function to calculate differences and confidence intervals
calculate_diff_ci <- function(model, term1, term2) {
  coefs <- coef(model)
  vcov_matrix <- vcov(model)
  
  estimate <- coefs[term1] - coefs[term2]
  std_error <- sqrt(vcov_matrix[term1, term1] + vcov_matrix[term2, term2] - 2 * vcov_matrix[term1, term2])
  
  t_value <- qt(0.975, df = model$df.residual)
  
  ci_lower <- estimate - t_value * std_error
  ci_upper <- estimate + t_value * std_error
  
  data.frame(term = paste(term1, "-", term2),
             estimate = estimate,
             conf.low = ci_lower,
             conf.high = ci_upper)
}

# Run regressions and calculate differences for all DVs
dvs <- c("dv_fav", "dv_coop", "dv_buy", "dv_weak")
labels <- c("Favorable View", "Support Cooperation", "Buy Products", "Russia Weak")
results <- lapply(dvs, run_regression, data = fig1data)
names(results) <- dvs

# Find the range of the estimates for consistent x-axis scaling
left_estimates <- unlist(lapply(results, function(result) {
  result$tidy_model$estimate
}))
right_estimates <- unlist(lapply(results, function(result) {
  c(calculate_diff_ci(result$model, "autshort", "autnone")$estimate,
    calculate_diff_ci(result$model, "demshort", "demnone")$estimate)
}))

x_range_left <- c(10, 90)  # Fixed range from 10 to 90 for left-hand side
x_range_right <- c(-5, 50)  # Fixed start at -5 for right-hand side

# Function to create plots for each DV
create_plots <- function(result, dv_name, adjust_label = FALSE) {
  levels_data <- result$tidy_model %>%
    mutate(group = case_when(
      term %in% c("autnone", "autshort") ~ "Non-Democracy",
      term %in% c("demnone", "demshort") ~ "Democracy"
    ))
  
  diff_aut <- calculate_diff_ci(result$model, "autshort", "autnone")
  diff_dem <- calculate_diff_ci(result$model, "demshort", "demnone")
  
  treatment_data <- bind_rows(diff_aut, diff_dem) %>%
    mutate(group = c("Non-Democracy", "Democracy"))
  
  vjust_values <- ifelse(adjust_label & levels_data$term %in% c("autnone", "demnone"), 1.5, -1)
  
  p1 <- ggplot(levels_data, aes(y = group, x = estimate, xmin = conf.low, xmax = conf.high, color = term, shape = term)) +
    geom_pointrange() +
    geom_text(aes(label = sprintf("%.1f", estimate), vjust = vjust_values), size = 5) +
    scale_color_manual(values = c("autnone" = "black", "autshort" = "#808080", "demnone" = "black", "demshort" = "#808080")) +  # Darker gray for hollow dots
    scale_shape_manual(values = c("autnone" = 16, "autshort" = 1, "demnone" = 16, "demshort" = 1)) +
    scale_x_continuous(breaks = seq(10, 90, by = 20)) +  # Set breaks every 20 units
    theme_bw(base_size = 15) +  # Smaller base size for titles
    labs(y = NULL, x = NULL) +
    theme(legend.position = "none", 
          plot.margin = margin(5, 5, 5, 5),
          strip.background = element_rect(fill = "gray"),
          strip.text = element_text(size = 10, face = "bold"),  # Smaller titles
          plot.title = element_text(hjust = 0.5),
          panel.grid.major.x = element_blank(),  # Remove vertical gridlines
          panel.grid.minor.x = element_blank(),
          axis.text.y = element_text(size = 15),
          axis.title.y = element_text(size = 15)) +
    ggtitle(paste0(dv_name, ": Levels")) +  # Remove extra space
    coord_cartesian(xlim = x_range_left)
  
  p2 <- ggplot(treatment_data, aes(y = group, x = estimate, xmin = conf.low, xmax = conf.high)) +
    geom_pointrange(shape = 1, color = "black") +  # Hollow round black dots
    geom_text(aes(label = sprintf("%.1f", estimate)), vjust = -1, size = 5) +
    theme_bw(base_size = 15) +  # Smaller base size for titles
    labs(y = NULL, x = NULL) +
    theme(plot.margin = margin(5, 5, 5, 5),
          strip.background = element_rect(fill = "gray"),
          strip.text = element_text(size = 10, face = "bold"),  # Smaller titles
          plot.title = element_text(hjust = 0.5),
          panel.grid.major.x = element_blank(),  # Remove vertical gridlines
          panel.grid.minor.x = element_blank(),
          axis.text.y = element_blank()) +
    ggtitle(paste0(dv_name, ": Effects")) +  # Remove extra space
    coord_cartesian(xlim = x_range_right)
  
  list(p1 = p1, p2 = p2)
}

# Create plots for each DV
plots <- mapply(create_plots, results, labels, adjust_label = c(FALSE, FALSE, FALSE, TRUE), SIMPLIFY = FALSE)

# Remove y-axis labels for the right-hand side plots
for (i in seq_along(plots)) {
  plots[[i]]$p2 <- plots[[i]]$p2 + theme(axis.title.y = element_blank())
}

# Combine all plots into a single figure
combined_plot <- plot_grid(
  plot_grid(plots[[1]]$p1, plots[[1]]$p2, ncol = 2, rel_widths = c(1, .75)),
  plot_grid(plots[[2]]$p1, plots[[2]]$p2, ncol = 2, rel_widths = c(1, .75)),
  plot_grid(plots[[3]]$p1, plots[[3]]$p2, ncol = 2, rel_widths = c(1, .75)),
  plot_grid(plots[[4]]$p1, plots[[4]]$p2, ncol = 2, rel_widths = c(1, .75)),
  ncol = 1,
  rel_heights = c(1, 1, 1, 1) * 0.15  # Further reduce vertical white space
)

# Save the combined plot in PNG format
png("figure1_main_analysis.png", width = 2400, height = 2400, res = 300)  # High resolution
print(combined_plot)
dev.off()

# Display the PNG in RStudio viewer
img <- readPNG("figure1_main_analysis.png")
grid::grid.raster(img)

# Save the combined plot in EMF format using devEMF
emf("figure1_main_analysis.emf", width = 12, height = 9)
print(combined_plot)
dev.off()

# Create revised and more compact legend with adjusted alignment
legend <- ggplot() + 
  geom_point(aes(x = 1, y = 2), color = "black", shape = 16, size = 5) + 
  geom_point(aes(x = 1, y = 1.5), color = "#808080", shape = 1, size = 5) + 
  geom_point(aes(x = 3, y = 2), color = "black", shape = 1, size = 5) + 
  theme_void() +
  annotate("text", x = 1.1, y = 2, label = "Level of DV with No Apology", hjust = 0, size = 5) +
  annotate("text", x = 1.1, y = 1.5, label = "Level of DV with Apology", hjust = 0, size = 5) +
  annotate("text", x = 3.1, y = 2, label = "Effect of Apology", hjust = 0, size = 5) +
  xlim(0.5, 4) + 
  ylim(1.3, 2.2) + 
  theme(legend.position = "none")

# Save the legend in EMF format using devEMF
emf("figure1_main_analysis_legend.emf", width = 12, height = 1)
print(legend)
dev.off()

########################### Figure 2 ###########################

fig2data <- read_dta("Study2_prepped.dta")

# Function to run regression and tidy results
run_regression <- function(dv, data) {
  model <- lm(as.formula(paste(dv, "~ autnoexpl + autnom + autshort + demnoexpl + demnom + demshort - 1")), data = fig2data)
  tidy_model <- tidy(model, conf.int = TRUE)
  tidy_model$term <- factor(tidy_model$term, levels = c("autnoexpl", "autnom", "autshort", "demnoexpl", "demnom", "demshort"))
  list(model = model, tidy_model = tidy_model)
}

# Function to calculate differences and confidence intervals
calculate_diff_ci <- function(model, term1, term2) {
  coefs <- coef(model)
  vcov_matrix <- vcov(model)
  
  estimate <- coefs[term1] - coefs[term2]
  std_error <- sqrt(vcov_matrix[term1, term1] + vcov_matrix[term2, term2] - 2 * vcov_matrix[term1, term2])
  
  t_value <- qt(0.975, df = model$df.residual)
  
  ci_lower <- estimate - t_value * std_error
  ci_upper <- estimate + t_value * std_error
  
  data.frame(term = paste(term1, "-", term2),
             estimate = estimate,
             conf.low = ci_lower,
             conf.high = ci_upper)
}

# Run regressions and calculate differences for all DVs
dvs <- c("dv_fav", "dv_coop", "dv_buy", "dv_weak")
labels <- c("Favorable View", "Support Cooperation", "Buy Products", "Russia Weak")
results <- lapply(dvs, run_regression, data = fig2data)
names(results) <- dvs

# Find the range of the estimates for consistent x-axis scaling
left_estimates <- unlist(lapply(results, function(result) {
  result$tidy_model$estimate
}))
right_estimates <- unlist(lapply(results, function(result) {
  c(
    calculate_diff_ci(result$model, "autshort", "autnoexpl")$estimate,
    calculate_diff_ci(result$model, "autshort", "autnom")$estimate,
    calculate_diff_ci(result$model, "demshort", "demnoexpl")$estimate,
    calculate_diff_ci(result$model, "demshort", "demnom")$estimate
  )
}))

x_range_left <- c(10, 80)  # Adjusted range from 10 to 80 for left-hand side
x_range_right <- c(0, 40)  # Adjusted range from 5 to 40 for right-hand side

# Function to create plots for each DV
create_plots <- function(result, dv_name, adjust_label = FALSE) {
  levels_data <- result$tidy_model %>%
    mutate(group = case_when(
      term %in% c("autnoexpl", "autnom", "autshort") ~ "Non-Democracy",
      term %in% c("demnoexpl", "demnom", "demshort") ~ "Democracy"
    ))
  
  diff_aut_noexpl <- calculate_diff_ci(result$model, "autshort", "autnoexpl")
  diff_aut_nom <- calculate_diff_ci(result$model, "autshort", "autnom")
  diff_dem_noexpl <- calculate_diff_ci(result$model, "demshort", "demnoexpl")
  diff_dem_nom <- calculate_diff_ci(result$model, "demshort", "demnom")
  
  treatment_data <- bind_rows(diff_aut_noexpl, diff_aut_nom, diff_dem_noexpl, diff_dem_nom) %>%
    mutate(group = rep(c("Non-Democracy", "Non-Democracy", "Democracy", "Democracy"), each = 1),
           shape = c(1, 22, 1, 22))  # Different shapes for different comparisons
  
  vjust_values <- case_when(
    levels_data$term %in% c("autnoexpl", "demnoexpl") ~ -1,  # Explicit no apology above the points
    levels_data$term %in% c("autnom", "demnom") ~ 1.5,       # No mention below the points
    TRUE ~ -1                                               # Apology treatments above the points
  )
  
  x_range <- if (dv_name == "Russia Weak") c(0, 30) else x_range_left  # Adjust x-range for dv_weak and show "0"
  
  p1 <- ggplot(levels_data, aes(y = group, x = estimate, xmin = conf.low, xmax = conf.high, color = term, shape = term)) +
    geom_pointrange() +
    geom_text(aes(label = sprintf("%.1f", estimate), vjust = vjust_values), size = 5) +
    scale_color_manual(values = c("autnoexpl" = "black", "autnom" = "#606060", "autshort" = "#808080", "demnoexpl" = "black", "demnom" = "#606060", "demshort" = "#808080")) +
    scale_shape_manual(values = c("autnoexpl" = 16, "autnom" = 17, "autshort" = 1, "demnoexpl" = 16, "demnom" = 17, "demshort" = 1)) +
    scale_x_continuous(breaks = seq(0, 80, by = 10)) +  # Set breaks every 10 units, starting at 0
    theme_bw(base_size = 15) +  # Smaller base size for titles
    labs(y = NULL, x = NULL) +
    theme(legend.position = "none",  # Remove legend from left panel
          plot.margin = margin(5, 5, 5, 5),
          strip.background = element_rect(fill = "gray"),
          strip.text = element_text(size = 10, face = "bold"),  # Smaller titles
          plot.title = element_text(hjust = 0.5),
          panel.grid.major.x = element_blank(),  # Remove vertical gridlines
          panel.grid.minor.x = element_blank(),
          axis.text.y = element_text(size = 15),
          axis.title.y = element_text(size = 15)) +
    ggtitle(paste0(dv_name, ": Levels")) +  # Remove extra space
    coord_cartesian(xlim = x_range)
  
  x_range_right_adjusted <- if (dv_name == "Russia Weak") c(-10, 10) else x_range_right  # Adjust x-range for dv_weak
  
  vjust_values_right <- c(-1, 1.5, -1, 1.5)  # Explicit no apology comparisons above the line, no mention comparisons below
  
  p2 <- ggplot(treatment_data, aes(y = group, x = estimate, xmin = conf.low, xmax = conf.high, shape = factor(shape))) +
    geom_pointrange(color = "black") +  # Hollow round black dots
    geom_text(aes(label = sprintf("%.1f", estimate), vjust = vjust_values_right), size = 5) +
    scale_shape_manual(values = c(1, 22)) +  # Different shapes for different comparisons
    theme_bw(base_size = 15) +  # Smaller base size for titles
    labs(y = NULL, x = NULL) +
    theme(legend.position = "none",  # Ensure no legends are shown in the right panel
          plot.margin = margin(5, 5, 5, 5),
          strip.background = element_rect(fill = "gray"),
          strip.text = element_text(size = 10, face = "bold"),  # Smaller titles
          plot.title = element_text(hjust = 0.5),
          panel.grid.major.x = element_blank(),  # Remove vertical gridlines
          panel.grid.minor.x = element_blank(),
          axis.text.y = element_blank()) +
    ggtitle(paste0(dv_name, ": Effects")) +  # Remove extra space
    coord_cartesian(xlim = x_range_right_adjusted)
  
  list(p1 = p1, p2 = p2)
}

# Create plots for each DV
plots <- mapply(create_plots, results, labels, adjust_label = c(FALSE, FALSE, FALSE, TRUE), SIMPLIFY = FALSE)

# Combine all plots into a single figure
combined_plot <- plot_grid(
  plot_grid(plots[[1]]$p1, plots[[1]]$p2, ncol = 2, rel_widths = c(1, .75)),
  plot_grid(plots[[2]]$p1, plots[[2]]$p2, ncol = 2, rel_widths = c(1, .75)),
  plot_grid(plots[[3]]$p1, plots[[3]]$p2, ncol = 2, rel_widths = c(1, .75)),
  plot_grid(plots[[4]]$p1, plots[[4]]$p2, ncol = 2, rel_widths = c(1, .75)),
  ncol = 1,
  rel_heights = c(1, 1, 1, 1) * 0.15  # Further reduce vertical white space
)

# Save the combined plot in PNG format
png("figure2_main_analysis.png", width = 2400, height = 2400, res = 300)  # High resolution
print(combined_plot)
dev.off()

# Display the PNG in RStudio viewer
img <- readPNG("figure2_main_analysis.png")
grid::grid.raster(img)

# Save the combined plot in EMF format using devEMF
emf("figure2_main_analysis.emf", width = 12, height = 9)
print(combined_plot)
dev.off()

# Create revised and more compact legend with adjusted alignment
legend <- ggplot() + 
  geom_point(aes(x = 1, y = 3), color = "black", shape = 16, size = 5) + 
  geom_point(aes(x = 1, y = 2.5), color = "#606060", shape = 17, size = 5) + 
  geom_point(aes(x = 1, y = 2), color = "#808080", shape = 1, size = 5) + 
  geom_point(aes(x = 3, y = 3), color = "black", shape = 1, size = 5) + 
  geom_point(aes(x = 3, y = 2.5), color = "black", shape = 22, size = 5) + 
  theme_void() +
  annotate("text", x = 1.1, y = 3, label = "Level of DV with Explicit No Apology", hjust = 0, size = 5) +
  annotate("text", x = 1.1, y = 2.5, label = "Level of DV with No Mention of Apology", hjust = 0, size = 5) +
  annotate("text", x = 1.1, y = 2, label = "Level of DV with Apology", hjust = 0, size = 5) +
  annotate("text", x = 3.1, y = 3, label = "Effect of Apology vs Explicit No Apology", hjust = 0, size = 5) +
  annotate("text", x = 3.1, y = 2.5, label = "Effect of Apology vs No Mention of Apology", hjust = 0, size = 5) +
  xlim(0.5, 5) + 
  ylim(1.8, 3.2) + 
  theme(legend.position = "none")

# Save the legend in EMF format using devEMF
emf("figure2_main_analysis_legend.emf", width = 12, height = 1)
print(legend)
dev.off()


####################################### Figure 3 ####################################

# Filter data to include only no apology and long apology without backlash
fig3data <- data_full %>% filter(apol %in% c(0, 3, 5))

# Function to run regression and tidy results
run_regression <- function(dv, data) {
  model <- lm(as.formula(paste(dv, "~ autnone + autreject + autaccept + demnone + demreject + demaccept - 1")), data = fig3data)
  tidy_model <- tidy(model, conf.int = TRUE)
  tidy_model$term <- factor(tidy_model$term, levels = c("autnone", "autreject", "autaccept", "demnone", "demreject", "demaccept"))
  list(model = model, tidy_model = tidy_model)
}

# Function to calculate differences and confidence intervals
calculate_diff_ci <- function(model, term1, term2) {
  coefs <- coef(model)
  vcov_matrix <- vcov(model)
  
  estimate <- coefs[term1] - coefs[term2]
  std_error <- sqrt(vcov_matrix[term1, term1] + vcov_matrix[term2, term2] - 2 * vcov_matrix[term1, term2])
  
  t_value <- qt(0.975, df = model$df.residual)
  
  ci_lower <- estimate - t_value * std_error
  ci_upper <- estimate + t_value * std_error
  
  data.frame(term = paste(term1, "-", term2),
             estimate = estimate,
             conf.low = ci_lower,
             conf.high = ci_upper)
}

# Run regressions and calculate differences for all DVs
dvs <- c("dv_fav", "dv_coop", "dv_buy", "dv_weak")
labels <- c("Favorable View", "Support Cooperation", "Buy Products", "Russia Weak")
results <- lapply(dvs, run_regression, data = fig3data)
names(results) <- dvs

# Find the range of the estimates for consistent x-axis scaling
left_estimates <- unlist(lapply(results, function(result) {
  result$tidy_model$estimate
}))
right_estimates <- unlist(lapply(results, function(result) {
  c(calculate_diff_ci(result$model, "autreject", "autnone")$estimate,
    calculate_diff_ci(result$model, "autaccept", "autnone")$estimate,
    calculate_diff_ci(result$model, "demreject", "demnone")$estimate,
    calculate_diff_ci(result$model, "demaccept", "demnone")$estimate)
}))

x_range_left <- c(10, 90)  # Fixed range from 10 to 90 for left-hand side
x_range_right <- c(-5, 58)  # Extend right-hand side range to 58

# Function to create plots for each DV
create_plots <- function(result, dv_name, adjust_label = FALSE) {
  levels_data <- result$tidy_model %>%
    mutate(group = case_when(
      term %in% c("autnone", "autreject", "autaccept") ~ "Non-Democracy",
      term %in% c("demnone", "demreject", "demaccept") ~ "Democracy"
    ))
  
  diff_aut_reject <- calculate_diff_ci(result$model, "autreject", "autnone")
  diff_aut_accept <- calculate_diff_ci(result$model, "autaccept", "autnone")
  diff_dem_reject <- calculate_diff_ci(result$model, "demreject", "demnone")
  diff_dem_accept <- calculate_diff_ci(result$model, "demaccept", "demnone")
  
  treatment_data <- bind_rows(diff_aut_reject, diff_aut_accept, diff_dem_reject, diff_dem_accept) %>%
    mutate(group = rep(c("Non-Democracy", "Non-Democracy", "Democracy", "Democracy"), each = 1),
           shape = c(1, 22, 1, 22))  # Different shapes for different comparisons
  
  vjust_values <- ifelse(adjust_label & levels_data$term %in% c("autnone", "demnone"), 1.5, -1)
  
  x_range <- if (dv_name == "Russia Weak") c(5, 25) else x_range_left  # Adjust x-range for dv_weak
  
  p1 <- ggplot(levels_data, aes(y = group, x = estimate, xmin = conf.low, xmax = conf.high, color = term, shape = term)) +
    geom_pointrange() +
    geom_text(aes(label = sprintf("%.1f", estimate), vjust = vjust_values), size = 5) +
    scale_color_manual(values = c("autnone" = "black", "autreject" = "#808080", "autaccept" = "#404040", "demnone" = "black", "demreject" = "#808080", "demaccept" = "#404040")) +  # Different shades of gray for different apology types
    scale_shape_manual(values = c("autnone" = 16, "autreject" = 1, "autaccept" = 2, "demnone" = 16, "demreject" = 1, "demaccept" = 2)) +
    scale_x_continuous(breaks = if (dv_name == "Russia Weak") seq(5, 25, by = 5) else seq(10, 90, by = 20)) +  # Set breaks every 5 units for dv_weak
    theme_bw(base_size = 15) +  # Smaller base size for titles
    labs(y = NULL, x = NULL) +
    theme(legend.position = "none", 
          plot.margin = margin(5, 5, 5, 5),
          strip.background = element_rect(fill = "gray"),
          strip.text = element_text(size = 10, face = "bold"),  # Smaller titles
          plot.title = element_text(hjust = 0.5),
          panel.grid.major.x = element_blank(),  # Remove vertical gridlines
          panel.grid.minor.x = element_blank(),
          axis.text.y = element_text(size = 15),
          axis.title.y = element_text(size = 15)) +
    ggtitle(paste0(dv_name, ": Levels")) +  # Remove extra space
    coord_cartesian(xlim = x_range)
  
  x_range_right_adjusted <- if (dv_name == "Russia Weak") c(-10, 10) else x_range_right  # Adjust x-range for dv_weak
  
  p2 <- ggplot(treatment_data, aes(y = group, x = estimate, xmin = conf.low, xmax = conf.high, shape = factor(shape))) +
    geom_pointrange(color = "black") +  # Hollow round black dots
    geom_text(aes(label = sprintf("%.1f", estimate)), vjust = -1, size = 5) +
    scale_shape_manual(values = c(1, 22)) +  # Different shapes for different comparisons
    theme_bw(base_size = 15) +  # Smaller base size for titles
    labs(y = NULL, x = NULL) +
    theme(plot.margin = margin(5, 5, 5, 5),
          strip.background = element_rect(fill = "gray"),
          strip.text = element_text(size = 10, face = "bold"),  # Smaller titles
          plot.title = element_text(hjust = 0.5),
          panel.grid.major.x = element_blank(),  # Remove vertical gridlines
          panel.grid.minor.x = element_blank(),
          axis.text.y = element_blank(),
          legend.position = "none") +  # Ensure no legends are shown in the right panel
    ggtitle(paste0(dv_name, ": Effects")) +  # Remove extra space
    coord_cartesian(xlim = x_range_right_adjusted)
  
  list(p1 = p1, p2 = p2)
}

# Create plots for each DV
plots <- mapply(create_plots, results, labels, adjust_label = c(FALSE, FALSE, FALSE, TRUE), SIMPLIFY = FALSE)

# Combine all plots into a single figure
combined_plot <- plot_grid(
  plot_grid(plots[[1]]$p1, plots[[1]]$p2, ncol = 2, rel_widths = c(1, .75)),
  plot_grid(plots[[2]]$p1, plots[[2]]$p2, ncol = 2, rel_widths = c(1, .75)),
  plot_grid(plots[[3]]$p1, plots[[3]]$p2, ncol = 2, rel_widths = c(1, .75)),
  plot_grid(plots[[4]]$p1, plots[[4]]$p2, ncol = 2, rel_widths = c(1, .75)),
  ncol = 1,
  rel_heights = c(1, 1, 1, 1) * 0.15  # Further reduce vertical white space
)

# Save the combined plot in PNG format
png("figure3_main_analysis.png", width = 2400, height = 2400, res = 300)  # High resolution
print(combined_plot)
dev.off()

# Display the PNG in RStudio viewer
img <- readPNG("figure3_main_analysis.png")
grid::grid.raster(img)

# Save the combined plot in EMF format using devEMF
emf("figure3_main_analysis.emf", width = 12, height = 9)
print(combined_plot)
dev.off()

# Create revised and more compact legend with adjusted alignment
legend <- ggplot() + 
  geom_point(aes(x = 1, y = 3), color = "black", shape = 16, size = 5) + 
  geom_point(aes(x = 1, y = 2.5), color = "#808080", shape = 1, size = 5) + 
  geom_point(aes(x = 1, y = 2), color = "#404040", shape = 2, size = 5) + 
  geom_point(aes(x = 3, y = 3), color = "black", shape = 1, size = 5) + 
  geom_point(aes(x = 3, y = 2.5), color = "black", shape = 22, size = 5) + 
  theme_void() +
  annotate("text", x = 1.1, y = 3, label = "Level of DV with No Apology", hjust = 0, size = 5) +
  annotate("text", x = 1.1, y = 2.5, label = "Level of DV with Rejected Apology", hjust = 0, size = 5) +
  annotate("text", x = 1.1, y = 2, label = "Level of DV with Accepted Apology", hjust = 0, size = 5) +
  annotate("text", x = 3.1, y = 3, label = "Effect of Rejected Apology Relative to No Apology", hjust = 0, size = 5) +
  annotate("text", x = 3.1, y = 2.5, label = "Effect of Accepted Apology Relative to No Apology", hjust = 0, size = 5) +
  xlim(0.5, 5) + 
  ylim(1.8, 3.2) + 
  theme(legend.position = "none")

# Save the legend in EMF format using devEMF
emf("figure3_main_analysis_legend.emf", width = 12, height = 1)
print(legend)
dev.off()


####################################### Figure 4 ####################################

# Reload the data
data <- data_full

# Revise function to run regression and tidy results
run_regression <- function(dv, data) {
  model <- lm(as.formula(paste(dv, "~ autnone + autbacklash + autnobacklash + demnone + dembacklash + demnobacklash - 1")), data = data)
  tidy_model <- tidy(model, conf.int = TRUE)
  tidy_model$term <- factor(tidy_model$term, levels = c("autnone", "autbacklash", "autnobacklash", "demnone", "dembacklash", "demnobacklash"))
  list(model = model, tidy_model = tidy_model)
}

# Function to calculate differences and confidence intervals
calculate_diff_ci <- function(model, term1, term2) {
  coefs <- coef(model)
  vcov_matrix <- vcov(model)
  
  estimate <- coefs[term1] - coefs[term2]
  std_error <- sqrt(vcov_matrix[term1, term1] + vcov_matrix[term2, term2] - 2 * vcov_matrix[term1, term2])
  
  t_value <- qt(0.975, df = model$df.residual)
  
  ci_lower <- estimate - t_value * std_error
  ci_upper <- estimate + t_value * std_error
  
  data.frame(term = paste(term1, "-", term2),
             estimate = estimate,
             conf.low = ci_lower,
             conf.high = ci_upper)
}

# Function to generate figures based on filter criteria and output prefix
generate_figure <- function(filter_criteria, output_prefix) {
  # Filter data based on criteria
  filtered_data <- data %>% filter(eval(parse(text = filter_criteria)))
  
  # Run regressions and calculate differences for all DVs except dv_weak
  dvs <- c("dv_fav", "dv_coop", "dv_buy")
  labels <- c("Favorable View", "Support Cooperation", "Buy Products")
  results <- lapply(dvs, run_regression, data = filtered_data)
  names(results) <- dvs
  
  # Find the range of the estimates for consistent x-axis scaling
  left_estimates <- unlist(lapply(results, function(result) {
    result$tidy_model$estimate
  }))
  right_estimates <- unlist(lapply(results, function(result) {
    c(
      calculate_diff_ci(result$model, "autbacklash", "autnone")$estimate,
      calculate_diff_ci(result$model, "autnobacklash", "autnone")$estimate,
      calculate_diff_ci(result$model, "dembacklash", "demnone")$estimate,
      calculate_diff_ci(result$model, "demnobacklash", "demnone")$estimate
    )
  }))
  
  x_range_left <- c(10, 90)  # Fixed range from 10 to 90 for left-hand side
  x_range_right <- c(-5, 58)  # Extend right-hand side range to 58
  
  # Function to create plots for each DV
  create_plots <- function(result, dv_name, adjust_label = FALSE) {
    levels_data <- result$tidy_model %>%
      mutate(group = case_when(
        term %in% c("autnone", "autbacklash", "autnobacklash") ~ "Non-Democracy",
        term %in% c("demnone", "dembacklash", "demnobacklash") ~ "Democracy"
      ))
    
    diff_aut_backlash <- calculate_diff_ci(result$model, "autbacklash", "autnone")
    diff_aut_nobacklash <- calculate_diff_ci(result$model, "autnobacklash", "autnone")
    diff_dem_backlash <- calculate_diff_ci(result$model, "dembacklash", "demnone")
    diff_dem_nobacklash <- calculate_diff_ci(result$model, "demnobacklash", "demnone")
    
    treatment_data <- bind_rows(diff_aut_backlash, diff_aut_nobacklash, diff_dem_backlash, diff_dem_nobacklash) %>%
      mutate(group = rep(c("Non-Democracy", "Non-Democracy", "Democracy", "Democracy"), each = 1),
             shape = c(1, 22, 1, 22))  # Different shapes for different comparisons
    
    vjust_values <- ifelse(adjust_label & levels_data$term %in% c("autnone", "demnone"), 1.5, -1)
    
    x_range <- x_range_left  # Adjust x-range
    
    p1 <- ggplot(levels_data, aes(y = group, x = estimate, xmin = conf.low, xmax = conf.high, color = term, shape = term)) +
      geom_pointrange() +
      geom_text(aes(label = sprintf("%.1f", estimate), vjust = vjust_values), size = 5) +
      scale_color_manual(values = c("autnone" = "black", "autbacklash" = "#808080", "autnobacklash" = "#404040", "demnone" = "black", "dembacklash" = "#808080", "demnobacklash" = "#404040")) +
      scale_shape_manual(values = c("autnone" = 16, "autbacklash" = 1, "autnobacklash" = 22, "demnone" = 16, "dembacklash" = 1, "demnobacklash" = 22)) +
      scale_x_continuous(breaks = seq(10, 90, by = 20)) +  # Set breaks every 20 units
      theme_bw(base_size = 15) +  # Smaller base size for titles
      labs(y = NULL, x = NULL) +
      theme(legend.position = "none", 
            plot.margin = margin(5, 5, 5, 5),
            strip.background = element_rect(fill = "gray"),
            strip.text = element_text(size = 10, face = "bold"),  # Smaller titles
            plot.title = element_text(hjust = 0.5),
            panel.grid.major.x = element_blank(),  # Remove vertical gridlines
            panel.grid.minor.x = element_blank(),
            axis.text.y = element_text(size = 15),
            axis.title.y = element_text(size = 15)) +
      ggtitle(paste0(dv_name, ": Levels")) +  # Remove extra space
      coord_cartesian(xlim = x_range)
    
    p2 <- ggplot(treatment_data, aes(y = group, x = estimate, xmin = conf.low, xmax = conf.high, shape = factor(shape))) +
      geom_pointrange(color = "black") +  # Hollow round black dots
      geom_text(aes(label = sprintf("%.1f", estimate)), vjust = -1, size = 5) +
      scale_shape_manual(values = c(1, 22)) +  # Different shapes for different comparisons
      theme_bw(base_size = 15) +  # Smaller base size for titles
      labs(y = NULL, x = NULL) +
      theme(plot.margin = margin(5, 5, 5, 5),
            strip.background = element_rect(fill = "gray"),
            strip.text = element_text(size = 10, face = "bold"),  # Smaller titles
            plot.title = element_text(hjust = 0.5),
            panel.grid.major.x = element_blank(),  # Remove vertical gridlines
            panel.grid.minor.x = element_blank(),
            axis.text.y = element_blank(),
            legend.position = "none") +  # Ensure no legends are shown in the right panel
      ggtitle(paste0(dv_name, ": Effects")) +  # Remove extra space
      coord_cartesian(xlim = x_range_right)
    
    list(p1 = p1, p2 = p2)
  }
  
  # Create plots for each DV
  plots <- mapply(create_plots, results, labels, adjust_label = c(FALSE, FALSE, FALSE), SIMPLIFY = FALSE)
  
  # Combine all plots into a single figure
  combined_plot <- plot_grid(
    plot_grid(plots[[1]]$p1, plots[[1]]$p2, ncol = 2, rel_widths = c(1, .75)),
    plot_grid(plots[[2]]$p1, plots[[2]]$p2, ncol = 2, rel_widths = c(1, .75)),
    plot_grid(plots[[3]]$p1, plots[[3]]$p2, ncol = 2, rel_widths = c(1, .75)),
    ncol = 1,
    rel_heights = c(1, 1, 1) * 0.15  # Further reduce vertical white space
  )
  
  # Save the combined plot in PNG format
  png(paste0(output_prefix, "_main_analysis.png"), width = 2400, height = 2400, res = 300)  # High resolution
  print(combined_plot)
  dev.off()
  
  # Display the PNG in RStudio viewer
  img <- readPNG(paste0(output_prefix, "_main_analysis.png"))
  grid::grid.raster(img)
  
  # Save the combined plot in EMF format using devEMF
  emf(paste0(output_prefix, "_main_analysis.emf"), width = 12, height = 9)
  print(combined_plot)
  dev.off()
  
  # Create revised and more compact legend with adjusted alignment
  legend <- ggplot() + 
    geom_point(aes(x = 1, y = 3), color = "black", shape = 16, size = 5) + 
    geom_point(aes(x = 1, y = 2.5), color = "#808080", shape = 1, size = 5) + 
    geom_point(aes(x = 1, y = 2), color = "#404040", shape = 22, size = 5) + 
    geom_point(aes(x = 2.6, y = 3), color = "black", shape = 1, size = 5) + 
    geom_point(aes(x = 2.6, y = 2.5), color = "black", shape = 22, size = 5) + 
    theme_void() +
    annotate("text", x = 1.1, y = 3, label = "Level of DV: No Apology", hjust = 0, size = 5) +
    annotate("text", x = 1.1, y = 2.5, label = "Level of DV: Apology with Backlash", hjust = 0, size = 5) +
    annotate("text", x = 1.1, y = 2, label = "Level of DV: Apology", hjust = 0, size = 5) +
    annotate("text", x = 2.7, y = 3, label = "Effect: Apology with Backlash vs No Apology", hjust = 0, size = 5) +
    annotate("text", x = 2.7, y = 2.5, label = "Effect: Apology with Backlsh vs No Apology", hjust = 0, size = 5) +
    xlim(0.5, 5) + 
    ylim(1.8, 3.2) + 
    theme(legend.position = "none")
  
  # Save the legend in EMF format using devEMF
  emf(paste0(output_prefix, "_main_analysis_legend.emf"), width = 12, height = 1)
  print(legend)
  dev.off()
}

# Generate figure 4a for accepted apologies
generate_figure("apol %in% c(0, 4, 5)", "figure4a")

# Generate figure 4b for rejected apologies
generate_figure("apol %in% c(0, 2, 3)", "figure4b")



